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^' ABSTRACT 

O . 

$—( , We study the linear m = 1 counter-rotating instability in a two-component, 

jrt ' nearly Keplerian disc. Our goal is to understand these slow modes in discs 

orbiting massive black holes in galactic nuclei. They are of interest not only be- 
cause they are of large spatial scale — and can hence dominate observations — 
but also because they can be growing modes that are readily excited by ac- 



> 

en , 

Cn . cretion events. Self-gravity being nonlocal, the eigenvalue problem results in 



a pair of coupled integral equations, which we derive for a two-component 
O ; softened gravity disc. We solve this integral eigenvalue problem numerically 

for various values of mass fraction in the counter- rotating component. The 



eigenvalues are in general complex, being real only in the absence of the 
counter-rotating component, or imaginary when both components have iden- 
tical surface density profiles. Our main results are as follows: (i) the pattern 
speed appears to be non negative, with the growth (or damping) rate be- 
ing larger for larger values of the pattern speed; (ii) for a given value of 
the pattern speed, the growth (or damping) rate increases as the mass in 
the counter-rotating component increases; (iii) the number of nodes of the 
eigenfunctions decreases with increasing pattern speed and growth rate. Ob- 
servations of lopsided brightness distributions would then be dominated by 
modes with the least number of nodes, which also possess the largest pattern 
speeds and growth rates. 

Key words: instabilities — stellar dynamics — celestial mechanics — galax- 
ies: nuclei 
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2 Gulati, Saini & Sridhar 
1 INTRODUCTION 

Most galactic nuclei are thought to host massive black holes and dense clusters of stars whose 



structure and kinernatics are correlated with glo bal galaxy properties (JGebhardt et al. 



Ferrarese fc Merritt 



2000 



Gebhardt et al 



1996 : 



2000l ). Such correlations raise questions of great 



interest related to galaxy formation and the growth of nuclear black holes. The nearby 
large spiral galaxy M31 has a n off-centered peak in a double-peaked brightness distribution 



aroun d its nuclear black hole ( jLight et al. 



1974 : 



Lauer et al.l 



1993 



1998 



Kormendy fc Bender 



19991 ). This lopsided brightness distribution could arise naturally if the a poapses o 



stellar orbits, orbiting the black hole, happened to be clustered together (JTremaine 



many 



19951 ). 



Since then, kinematic and dynamical rnodels of such an eccentric disc have been con 



structed by several authors 



2002 



Peiris fc Tremaine 



Bacon et al. 



2001 



Salow &: Statler 



2001 



Sambhus fc Sridhai 



Sambhus fc Sridhar 



20031 ). Of particular interest to this work is the model ofy 
(J2002l ). which included a few per cent of stars on retrograde (i.e. counter-rotating) orbits. 
They proposed that these stars could have been accreted to the centre of M31 in the form 
of a globular clu ster tha t spira led in due to dynamical friction. This proposal was motivated 



by the work of 



Toumal (120021 ) . who demonstrated that a Keplerian axisymmetric disc is 



susceptible to a linear lopsided instability in the m = 1 mode, even when a small fraction 
of the disc mass is in retrograde motion. 



Toumal (|2002| ) considered the linearized secular dynamics of particles orbiting a point 
mass, wherein particle orbits may be thought of as slowly deforming elliptical rings of small 
eccentricities. The m = 1 counter-rotating instability was studied analytically for a two-ring 
system, and numerically fo r a many-ring system. T he corresponding problem for continuous 



discs was then studied by 



Sridhar fc Saini 



(I2OIOI ). who proposed a simple model with dy- 
namics that could be studied largely analytically in the Wentzel-Kramers-Brillouin (WKB) 
approximation. Their model consisted of a two-component softe ned gravity d isc, orbiting 
a massive central black hole. Softened gravity was introduced by iMillen (jl97ll ) to simplify 
the analysis of the dynamics of stellar systems. In this form of interaction, the Newtonian 



1/d gravitational potential is replaced by l/y/d'^ + 6^, where 6 > is called the softening 
length. In the context of waves in discs, it is well kn own that the softening len gth mimics 
the epicyclic radius of stars on nearly circular orbits ( iBinney &: Tremainell2008l ). Therefore, 
a disc composed of cold coUisionless matter, interacting via softened gravity, provides a 
surrogate for a hot coUisionless disc. 
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Sridhar fc Saini 



(l2010f ) used a short-wavelength (WKB) approximation, derived analyt- 
ical expressions for the dispersion relation and showed that the frequency u is smaller than 
the Keplerian orbital frequency by a factor proportional to the small quantity e = Mct/M 
(which is the ratio of the disc mass to mass of the central object); in other words, the modes 
are slow. The WKB dispersion relation was used to argue that equal mass counter-rotating 
discs with the same surface density profiles (i.e. when there is not net rotation) could have 
unstable modes. They also argued tha t, for an arbitrary mass ratio, the discs must be unre- 



alistically hot to avoid an instability. ISridhar fc Saini I ( 120101 ) then used Bohr-Sommerfeld 



quantization to construct global modes, within the WKB approximation. A matter of con- 
cern is that the wavelengths of the modes could be of order the scale length of the discs; 
the modes being large-sca le it is possible t h at th e WKB approximation itself is invalid. 



Another limitation is that 



Sridhar fc Saini 



( I2OIOI ) could construct (WKB) global modes 



only for the case of equal mass discs. Therefore it is necessary to address the full eigen- 
value problem to understand the systematic behaviour of eigenvalues and eigenf unctions. 
To this end, we formulate the eigenvalue problem for the linear, slow, m = 1 modes in a 
two-component, softened gravity, counter-rotating disc. Due to the long-range nature of 
gravitational interactions, we have to deal with a pair of coupled integral equations defining 
the eigenvalue problem. We draw some general conclusions and then proceed to solve the 
equations numerically for eigenvalues and eigenfunctions. 

In § 2 we introduce the unperturbed two-component nearly Keplerian disc, define the 
apse precession rates, discuss the potential theory for softened gravity, and derive the cou- 
pled, linear integral equations that determine the eigenvalue problems for slow m = 1 modes. 
A derivation of the relationship between the softened Laplace coefficients (used in the po- 
tential theory of § 2 and § 3) and the usual (unsoftened) Laplace coefficients is given in 
the Appendix. We specialize to discs with similar surface density profiles in § 3, when the 
two coupled equations can be cast as a single integral equation in a new mode variable; this 
allows us to draw some general conclusions about the eigenvalue problem. We also discuss 
in detail the numerical method to be employed. Our results are presented in § 4, where 
the properties of the stable, unstable and overstable modes are discussed. Conclusions are 
offered in § 5, where we seek to provide a global perspective on the correlations that occur 
between the pattern speeds, growth rates and eigenfunctions, as well as the variations of 
these quantities on the mass fraction in retrograde orbits. 
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2 FORMULATION OF THE LINEAR EIGENVALUE PROBLEM 

We consider linear non-axisymmetric perturbations in two counter-rotating discs orbiting a 
central point mass M. The discs are assumed to be coplanar and consist of cold collisionless 
particles which attract each other through softened gravity. However, the central mass and 
the disc particles interact via the usual (unsoftened) Newtonian gravity. Softened gravity is 
known to mimic the effects of velocity dispersion, so our discs are surrogates for hot stellar 
discs. We assume that the total mass in the discs, M^, is small in comparison to the central 
mass. Since e = Md/M ^ 1, the dynamics is dominated by the Keplerian attraction of the 
central mass, and the self-gravity of the discs is a small perturbation which enables slow 
modes. Below we formulate the linear eigenvalue problem of slow modes. 



2.1 Unperturbed discs 

We use polar-coordinates r = (-R , 0) in the plane of the discs, with the origin at the loca- 
tion of the central mass. Throughout this paper the superscripts '-I-' and ' — ' refer to the 
prograde and retrograde components, respectively. The unperturbed discs are assumed to 
be axisymmetric with surface densities S^ (R). The disc particles orbit in circles with veloci- 
ties, v^ = ± RQ{R)e^, where f2(i?) > is the angular speed determined by the unperturbed 
gravitational potential. 



^R) = -™ + UR) ■ (1) 

The first term on the right side is the Keplerian potential due to the central mass, and $d(-R) 
is the softened gravitational potential due to the combined self-gravities of both the discs: 



$,(r) = -G I Sl^2±S£2dV, (2) 

Ir — r'l + b"^ 



where b is the Miller softening length; the potential $d(-R) is 0(e) compared to GM/R. Test 
particles for nearly circular prograde orbits have azimuthal and radial frequencies, Q and k, 
given by 
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n\R) 



K^R) 



GM 1 d^rf 

+ 



R^ 

GM 



R di?' 
3d$rf 



d=^$. 



(3) 



(4) 



i?3 ' i? di? di?2 ■ 

The line of apsides of a nearly circular eccentric particle orbit of angular frequency ±Q{R), 
subjected only to gravity, precesses at a rate given by ±zo{R), where 



w{R) = n{R) - k{R) 



2 d 
RdR 



<l>d(i?) + Oie' 



(5) 



2n{R) \RdR ' di?2 
The cancellation of the 0(1) term, {GM/ R^), which is common to both Vl? and h? makes 

w ~ 0{e). This is the special feature of nearly Keplerian discs which is responsible for the 

existence of (slow) modes whose eigenfrequencies are ~ 0{e) when compared with orbital 

frequencies. 



2.2 Perturbed discs 

Let vj(r,t) = ■u^(r,t)e/? + t>^(r,t)e0 and S^(r,t) be infinitesimal perturbations to the 
velocity fields and surface densities of the ± discs, respectively. These satisfy the following 
linearized Euler and continuity equations: 



-^ + (v,^ ■ V)v,± + (v± ■ V)v| 



-v$„, 



^ + V ■ (S,^v,^ + S^,^ 



0, 



(6) 



(7) 



where $a(i',^) is the perturbing potential. Fourier analyzing the perturbations in t and 0, 
we seek solutions of the form, Xa(r, t) = J2m -^Ti^) exp[i(m0 — ut)] . Then 



u 



mi: 



D± 



d 
{±mil — w)— — ± 
dR 



,m± 



D± 



K^ d m 



i(±mf]-a;)S;^ 
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m± 



1 d 

RdR 



it„,»n±\ 



iR^,< 



2mVL 
R 


<T)™ 




(8) 


-u) 


^m 




(9) 


ir 

+ 1 


^2± m± _ 
T d a 


= 0, 


(10) 
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where 



D, 



± 



K^ — {±.mVt — oj) 



(11) 



The above equations determine m" 



r7i± „,m± 



and S^^ in terms of the perturbing potential $^ ; 



this would be the solution were $^ due to an external source. 

We are interested in modes for which $JJ^ arises from self gravity. In this case it depends on 
the total perturbed surface density, [SJ^^(i?) + S™~(_R)] . Manipulating the Poisson integral 
given in Eq. [21 we obtain 



K{R) 



R'dR'PUR,R') [sr(^') + sr (^')] 



(12) 



where the kernel 



PUR,R') =-^B[-\a,P) + "^^ 



(5m,l + Sm,-l) ■ 



(13) 



The second term on the right side is the indirect term arising from the fact that our coordi- 
nate system can be a non-inertial frame, because its origin is located on the central mass. The 
first term is the direct term coming from the perturbed self-gravity. Here i?< = min(i?, R') , 
i?> = max(_R, R') , a = R</ Ry and /3 = 6/-R> . The functions. 



i?f)(a,/3) =1 Tde- 
^ Jo (1 



cos m6 



2apsl+ a2j^2)./2 ' (^^) 

are "softened Laplace coefficients", introduced in iToumal (l2002l ). They can be expressed 
in terms of the usual (unsoftened) Laplace coefficients, as shown in Appendix A. We note 
that the unperturbed disc potential ^^ can be obtained from the unperturbed disc density, 
S+(i?) + S^ (i?) , by using Eq. ([12]) with m = . 



2.3 Slow m = 1 modes 

Modes with azimuthal wavenumber m = ±1 are slow in the sense that their eigenfrequencies, 
u, are smaller than the orbital frequency, Q, by a factor ^ 0(e) . Without loss of generality 
we may choose m = 1. In the slow mode approximation (JTremaind [200ll ) . we use the fact 
that fi > tj in Eqs. ([8])— ([U]), and write 
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u 



i± 



^ Df Y'dR 



' d 2fi 



,i± 



± 



Df 



VL d Vt 

2"Zr ^ 'r 



$1 



$i 



where 



±-^^K^ + ^^{RT^Wa^) + ^^^A^ 



Df = ±2n{uj T ro). 



0, 



(15) 
(16) 
(17) 



(18) 



Eqs. flT^ and flT^ imply the following relations between the perturbed velocity amplitudes: 



u 



i± 






2i^^ , 






m 



We use Eqs. ( !T9|) in the continuity equation (IT7|) to eliminate m^^ and write, 



± OEi^ 



2 d 



Combining Eqs. (112]), ([HD and (118])— (EO]) we obtain 



(R'^'^K^) ■ 



(20) 



[C^ T MR)]v'a^iR) 



°° di?'i?'i/2 r 5 



2R^n{R') \dR 
d 



[R^Pi{R,R' 



dR 



, [i?'^/^Sl(i?>^(i?0 - i2'V2s^-(i?>i-(i?')] . (21) 



We rewrite this by defining 



zHR) 



n{R) 



1/2 



,1± 



(22) 



use the fact that D.{R) oc R ^^^ for a Keplerian flow, and integrate by parts to obtain, 
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/■OO 1 D/ 

[uj T vj{R)]z^{R) = - j ^2TiR,R') 



+ / —2J^iR,R') 



j:j{R')j:^{Ry 
n{R')n{R) 

n{R')n{R) 



1/2 



z+{R 



1/2 



z-iR'), (23) 



where 



TiR, R!) 



1 



1 d 



l + l^^\Pi{R,B!) 



(24) 



2d\nR'J \' ' 2d\nR^ 
It is convenient to write S^(i?) = ri{R)T,^{R) and S^(i?) = (1 — ri{R))T,^{R) , where r]{R) 
is the local mass fraction in the unperturbed counter- rotating component; by definition, 
^ ?7(-R) ^ 1- Then, Eq. fl23|) can be recast as 



uz+{R) = +UJZ^{R) + 



dR 



R 



r[(i - ^(^'))(i - viR))]'^'}C{R, R')z^{R') 



dR 



R 



-[V{R'){1-V{R))Y/')C{R,R')Z~{R') 



UZ (i?) = -UJZ (R) + 



dR 



R 



-[{l-v{R')MR)]'^'}C{R,R')z^R') 



dR 



R 



yUR'UR)Y"lC{R,R')z-{R!) 



where the kernel 



(25) 



JC{R,R') 



1/2 



- 2 



2nG 



^AR')^<iiR) 
_ n{R')n{R) 

- ^AR')^ARy "" 
n{R')n{R) 

^AR')^AR) 
n{R')n{R) 



J^{R, R') 



i.i- ' 



2d\nR' 



l + l-J-^\PiiR,R') 



2(9 In i? 



.1/2 



?(1). 



' + 2 9h^J V + 2dh[R)^r^r~ • ^''^ 

Therefore the kernel /C(i?, i?') is a real symmetric function of R and i?'|j 

Using Eqs. ( IT9|) and ( l22l) . we can relate the eigenfunctions ;z'''(-R) and z~{R) to each 
other: 



^ The contribution from the indirect term in Pi{R, R') vanishes, because (2 + d/d In R') R' ^ = . 
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^{l-r]{R))[u-w{R)]z+{R) = ^/^[u + w{R)]z-{R) . (27) 

This relation can, in principle, be used to eliminate one of z'^{R) or z^{R) from the coupled 
Eqs. (|25l) . in which case the eigenvalue problem can be formulated in terms of a single func- 
tion (which can be either z'^{R) or z~{R)). However, such a procedure results in a further 
complication: the eigenvalue, u, will then occur inside the R' integral in the combination, 
{u ±w)/{uj ^ w), and this makes further analysis difficult. Eqs. (!25|) are symmetric under 
the (simultaneous) transformations, {'+' ,ri{R) ,00} — )■ {'— ' , [1 — //(i?)] ,—uj}, which in- 
terchange the meanings of the terms prograde and retrograde. It seems difficult to obtain 
general results when E^(i?) and S^(i?) have different functional forms. Below we consider 
the case when the mass fraction, rj, is a constant; i.e. when both Sj'(i?) and S^(i?) have the 
same radial profile. 



3 THE EIGENVALUE PROBLEM FOR CONSTANT r] DISCS 

When the counter-rotating discs have the same unperturbed surface density profiles, i.e. 
S^(i?) oc S^(i?) , some general results can be obtained. This case corresponds to the choice 
rj = constant, so that S^(_R) = rj'E^lR) and Sj"(i?) = {l — rj)'E^{R) . Then the eigenfunctions 
z~^{R) and z~{R) are related to each other by, 

[u-w{R)]^z+{R) = [uj + w{R)]^/{l-ri)z-{R). (28) 

Let us define a new function, Z{R), which is a linear combination of z^{R) and z~{R) : 



Z{R) = ^l-r]z+{R) - ^z-{R). (29) 

Then equations ( 12^ can be manipulated to derive a closed equation for Z{R) : 



w^ — w"^ 



/CO A T3I 
— ICiR,R')ZiR'), (30) 



[l-2r])u + ta_ 

We note that, in this integral eigenvalue problem for the single unknown function Z{R) , the 
(as yet undetermined) eigenvalue u occurs outside the integral. Once the problem is solved 
and Z{R) has been determined, we can use Eq. f l28|) and f l29|) to recover z^{R) : 
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(1 — 277)0; + tJ7 (1 — 2?7)CJ + tJ7 

Some general conclusions can be drawn: 

(i) In Eq. ( |30ll . the kernel )C{R,R') is real symmetric. Therefore, the eigenvalues, u, are 
either real or come in complex conjugate pairs. 

(ii) When n = 0, the counter-rotating component is absent, which is the case studied by 



Tremaind (1200 if ): then the left side of Eq. (!30|) becomes {uj — zo) Z . Since the kernel 1C{R, R') 
is real symmetric, the eigenvalues uj are real, so the slow modes are stable and oscillatory in 
time. Then the eigenfunctions, Z{R) may be taken to be real. Therefore z^{R) = Z{R) is a 
real function, and z~(R) = . g 

(iii) When rj = 1/2, there is equal mass in the counter-rotating component, and the 
surface densities of the ± discs are identical to each other. This case may also be thought 
of as one in which there is no net rotation at any radius. Then Eq. (|30|) becomes 



w-\R){u^-w\R))Z{R) = / ^}C{R,R')Z{R'), (32) 

Jo R 

Since the kernel /C(i?, R') is real symmetric, w^ must be real. There are two cases to consider, 
when the eigenvalues, w, are either real or purely imaginary. 

• When UJ is real, the slow modes are stable and oscillatory in time. The eigenfunctions 
z'^{R) can be taken to be real functions. 

• When u is imaginary, the eigenvalues come in pairs that are complex conjugates of each 
other, corresponding to non-oscillatory growing/damped modes. Let us set r] = 1/2 and 
u = i'j (where 7 is real) in Eq. (ISTl) : 

The function Z{R), which is a solution of Eq. ( l32|) . can be taken to be a real function 
multiplied by an arbitrary complex constant. It is useful to note two special cases: (i) when 
Z{R) is purely imaginary, then z~^{R) and z~{R) are complex conjugates of each other; (ii) 
when Z{R) is real, z~^{R) is equal to minus one times the complex conjugate of z~{R). 

To make progress for other values of r^, it seems necessary to address the eigenvalue 
problem numerically; the rest of this paper is devoted to this. 

■^ When 'rj = 1, the eigenvalues, oj, are again real, with z^ (R) = Z{R) a real function, and z'^{R) = . 
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3.1 Application to Kuzmin discs 

For numerical explorations of the eigenvalue problem, we consider the case when both the 
unperturbed ± discs are Kuzmin discs, with similar surface density profiles: S^(_R) = rjT,^{R) 
and S+(i?) = (1 - r])T.^{R) , where 



^AR) 



aM, 



(34) 



27r(i?2 + a2)3/2 

is the total surface density, Md is the total disc mass and a is the disc scale length. Kuzmin 
discs, being centrally concentrated, are re asonable candidates for unperturbed dis cs. More- 



over, earlier investigations of slow modes ( iTremaine 



2001 



Sridhar fc Saini 



2010h have ex- 



plored modes in Kuzmin discs, so we find this choice useful for comparisons with ear- 
lier work. The characteristic values of orbital frequency and surface density are given by 



VL* = ^GM ja^^ and S^ = Md/a^, respectively. The coupled Eqs. ( I25l) can be cast in a di- 
mensionless form in terms of these physical scales. The net effect is to rescale the eigenvalue 
u to a, where 



a 



Q*a 



u 



(35) 



In the following section, all quantities are to be taken as dimensionless; however, with some 
abuse of notation, we shall continue to use the same symbols for them. 



3.2 Numerical method 



Our method is broadly similar to iTremaind (1200 ll ). In order to calculate eigenvalues and 



eigenf unctions numerically, we approximate the integrals by a discrete sum using an Ap- 
point quadrature rule. The presence of the term dR/i? in the integrals suggests that a 
natural choice of variables is m = log(-R) and v = log(-R'), where, as mentioned above, R 
stands for the dimensionless length R/a. The use of a logarithmic scale is numerically more 
efficient, because it induces spacing in the coordinate space that increases with the radius. 
This handles naturally a certain expected behaviour of the eigenf unctions: since the surface 
density in a Kuzmin disc is a rapidly decreasing function of the radial distance, we expect 
the eigenfunctions to also decrease rapidly with increasing radius. Therefore, discretization 
of the coupled equations fl25|) follows the schema: 
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(36) 



where Wj's are suitably chosen weights. Then the discretized equations can be written as a 
matrix eigenvalue problem: 



where 



AC = ^C: 



(37) 



7])wjlCij + Wj5i_ 



■\/ri{l -vi)wjlCi 



^/riiY^^WjlCi 



-rjWjK,, 



«j 



Wjbij 



and 



V 



z- 



z, 



(38) 



The 2N X 2N matrix A has been represented above in a 2 x 2 block form, where each of 
the 4 blocks is a A^ x A^ matrix, with row and column indices i and j. Note that 5ij is the 
Kronecker delta symbol, and no summation is implied over the repeated j indices. Thus we 
have an eigenvalue problem for eigenvalues a, and eigenvectors given by the 2N dimensional 
column vector C • The use of unequal weights destroys the natural symmetry o f the kernel, 



but t 



lis is readily restored through a simple transformation given in § 8.1 of 



Press et al. 



(Il992[ ). The grid for our numerical calculations covers the range —7 ^ logi? ^ 5, which is 
divided into A^ = 4000 points; large r values of N give similar results. 



We note some differences with 



Tremaind (l200l[ ) concerning details of the numerical 



Tremaine 



metho d and assumptions. The major difference is in the treatment of softening: In 
teOOll ). a dimensionless softening parameter /3 = h/R was introduced, and the eigenvalue 
problem for slow modes was solved by holding the parameter /3 constant. This renders the 
physical softening length, 6, effectively dependent on radius, making it larger at larger radii, 
thereby not corresponding to any simple force law between two disc particles located at 
different radii. We have preferred to keep h constant, so that the force law between two disc 
parti cles is through th e usual Miller prescription. Other minor differences in treatment are: 



1 m 



Tremaind ( 120011 ). the disc interior to an inner cut-off radius was assumed to be frozen. 



In contras t we u se a straightforward inner cut-off radius of 10 ^ , as mentioned above; (ii) 
Tremaind ( 120011 ) uses a uniform grid in logi?, with four-point quadrature in the intervals 
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a = -0.2370 
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Figure 1. Slow g-modes in a single, prograde (77 = 0) Kuzniin disc with A = 0.1 and b = 10~ . The panels are labeled by the 
scaled eigenvalue cr. 



between consecutive grid points; we also use a uniform grid in logi? but instead employ a 
single A^-point quadrature for integration. 

Were we dealing with unsoftened gravity (i.e. the case when 6 = 0), the diagonal elements 
of the kernel would be singular. Hence, when the softening parameter b is much smaller than 
the grid size, accuracy is seriously compromised by round-off errors. Typically, the usable 
lower limit for b is ~ 10 ~^. 



4 NUMERICAL RESULTS 



We obtained the e igenvalues and eigen: 



'unctions of equation ( 1371) using the linear algebra 
package LAPACK (jAnderson et al.lll999h . We now present the results of our calculations for 
specific values of rj. As noted earlier, interchanging the meaning of prograde and retrograde 
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Figure 2. Slow p-modes in single, prograde (77 = 0) Kuzniin disc with no external source (i.e A = 1). The panels are labeled 
by the scaled eigenvalue a, and the softening parameter, b (which has been scaled with respect to a, the disc scale length). 



orbits leave the results invariant under the transformation (^7,0;) — )■ (1 — r/, — w); therefore, 
we present results below only for ^ r/ ^ 1/2 . 



4.1 No counter-rotation: 77 = 



es rotate i n the prograde sense. The eigen- 



Tremaind (1200 ll ). who also showed that the 



We are dealing with a single disc whose partic 

value problem for this case was studied first by 

eigenvalues are real; in other words, the disc supports stable slow modes. We consider this 

case first to benchmark our numerical method as well as assess the differences in results 

that may arise due to the mann er in which softe ning is treated. To facilitate comparison 



we use the same nomenclature as 



Tremaind (J200lh . Briefly, modes corresponding to positive 



and negative eigenvalues are referred to as "p-modes" and "g-modes" , respectively; we also 
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Figure 3. Growth rate versus number of nodes for r) = 0.5, for two values of softening, b = 0.1 and b = 0.05. 

introduce a parameter A = (1 + /)^^, where / is a constant that mimics additional precession 
due to an external source of the form tJ7e(_R) = fwdi^R)] we define eccentricity e^ through 



ex 



GM 



-1/2 



(39) 



and use the normalization, 



dR 



eliR) 



1 



(40) 



Our results corresponding to g-modes, for A = 0.1 and /3 = 10~^, are presented in 
Fig. [H where we plot modes with three or fewer nodes and give their eigenvalues. Results 
for p-modes for A = 1 (no external source) and various values of softening param eter b, are 



displayed in Fig. [2l These figures are to be compared with Fig. 3 and Fig. 6 of 
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Figure 4. Eigcnfunctions Z{R) are plotted as a function of the radial coordinate R, for rj = 0.5 . The panels are labeled by 
the values of the growth rate, T, and softening, h. 





Figure 5. Gray-scale plots of surface density perturbations, E^ (R,<j>,t) at time t = for the parameter values, r] = 0.5 and 
fe = 0.1 , and r = 0.0230. White/black correspond to the maximum positive/negative values of the perturbations. 



© 0000 RAS, MNRAS 000, 000-000 



Unstable m = 1 modes of counter-rotating Keplerian discs 17 



m 



(120011): the e j genf u nctions are of broadly similar form, but the eigenvalues differ from those 



Tremaind (120011 ) by upto ~ 30% . 



4.2 Equal counter-rotation (or no net rotation): rj = 1/2 

This case was studied by 



Sridhar fc Saini 



(]2010f ). who derived the following local or "WKB" 



dispersion relation: 



U) 



zu w 



n{R) 



k\ exp{—\k\b) 



(41) 



From this expression they concluded: if zu happens to be positive then u is real (and the 
disc is stable), but zu is negativ e for most continuous d iscs which implies that u can be 



either real or purely imaginary. 



Sridhar fc Saini 



( I2OIOI ) also studied global modes using 



Bohr-Sommerfeld quantization, which will be discussed later in this section. 

We have proved in the last section that the eigenvalues are either real (stable oscillatory 
modes), or purely imaginary (non-oscillatory growing and damped modes). Here we focus 
on the growing (unstable) modes; we define the growth rate of perturbations as. 



in order to facilitate comparisons with 



b\(T\ 



(42) 



Sridhar fc Saini 



mm. In Fig. IS l we p lot F for b 



Sridhar &: Saini 



torn found two 



equal to 0.1 and 0.05, versus the number of nodes of Z(R). 
separate branches in the spectrum, corresponding to long and short wavelengths (for each 
value of b). Comparing our Fig. |3]with their Figs. 3 & 4, we see that our results are more 
consistent with their short-wavelength branch than with their long-wavelength branch. This 
disagreement is probably because the long-wavelength branch corresponds to kR ~ 1, where 
WKB approximation breaks down. Moreover, the agreement between our results and their 
short-wavelength branch holds only in a broad sense, because there are differences in the 
numer ical values of the eigenvalues. We trace this difference to the fact that ISridhar fc Saini 
( 120101 ) used an analytical result for the precession frequency corresponding to unsoftened 
gravity, whereas we have consistently used softened gravity for all gravitational interactions 
between disc particles. This probably also results in another difference between our results: 



according to 



Sridhar fc Saini 



( I2OIOI ). 10"'^ < F < 10~^; however, as can be seen from our 
Fig. [31 we obtain values of F both inside and outside this range. Changing the value of b 
causes a horizontal shift in the spectrum, which is consistent with their results. 
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Figure 6. Distribution of eigenvalues in the complex a plane, for r] = 0.1,0.25 and 0.4. Panels labeled a give an overview, 
whereas the panels labeled b provide an close-up view near the origin. 



In Fig m we plot a few of these eigenfunctions as a function of R for both values of 
b equal to 0.1 and 0.05. Note that, from the discussion in the previous section, Z{R) can 
always be chosen to be a real function of R. The smallest number of nodes corresponds 
to the largest value of the growth factor. The eigenfunctions with the fewest nodes have 
significant amplitudes in a small range of radii around -R ~ 1, and this range increases with 
the number of nodes (and correspondingly, the growth rate decreases). Figgis a gray-scale 
plot of the surface density perturbations in the ± discs, S^(/2, 0, t = 0). Note the relative 
phase shift between the it perturbations. For other value of b shown in Fig H] we get similar 
patterns. 

© 0000 RAS, MNRAS 000, 000-000 



T 1 I I I I III 1 1 I I I INI 1 1 I I I INI 1 1 I I I I 1 11 




2 F 



1 

N 



■1 - 



-2 t 



n — I II Mill 



1 - 



N3 



-1 



-2 



Unstable m = 1 modes of counter-rotating Keplerian discs 
7) = 0.1 



19 



£7^=0.2501 
CTj=0.1754 



n — I I I I Mil 1 — I I I I Mil 1 — I I I 1 1 III 1 — I I I 1 1 1 11 



CTr=0.2221 
CTj=0.1524 ^ 




I I I MM 1 1 I I I I 1 11 



CTr=0.1574 
CTj=0.2635 



-q = 0.25 



n — I I 1 1 MM 




I I I 1 1 III 1 — I I I I II 11 



(Tr=0.1408 
aj=0.2332 



n 1 I I I I III 1 1 I I I MM 1 1 I I I MM 1 1 I I I I 1 11 



77 = 0.4 



CTp=0.0631 
CTj=0.3011 



r 



n — I I 1 1 MM 




— I I 1 1 1 ii| 1 — I I 1 1 II ij. 

(7(^=0.0565 
(Tj=0.2670 ^ 



qL 



0.1 



10 100 0.1 

Radial coordinate R 



10 100 



Figure 7. Real parts of the "most unstable" eigenfunctions Z{R), plotted as a function of the radial coordinate, R for b = 0.1 
and for r] = 0.1, 0.25 and 0.4. Panels are labeled by the real and imaginary parts of the eigenvalues. 



4.3 Other values of r] 

We present results for values of rj other t han and 1/2. T 



only because they were not explored by 



Sridhar fc Saini 



l ese ar e particularly interesting, not 
(J2OIOD, but because the eigenvalues 



can be truly complex, corresponding to growing and damped modes which precess with 
steady pattern speeds. We write the eigenvalues as a = aR + io"/. In Fig|6], we display the 
eigenvalues in the complex a-plane, for softening parameter 6 = 0.1 and for rj equal to 0.1, 
0.25 and 0.4. Panels on the left, labeled (a), provide an overview, whereas the panels on the 

ues n ear the origin 



right, labeled (6), provide a close-up view of the distribution o f eigenva 

of the complex a-plane; this distribution is similar to Fig. 3 of iToumal ( 20021 ). We are able 
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Figure 8. Gray-scale plots of surface density perturbations, E^ (_R, <f>, t) at time t = for the parameter values, r] = 0.25 and 
b = 0.1 , and cr = 0.1408 + iO.2332. White/black correspond to the maximum positive/negative values of the perturbations. 
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Figure 9. Real parts of two pairs of eigenfunctions Z{R) (from two arms of a branch), plotted as a function of the radial 
coordinate R, for b = 0.1 and i) = 0.25. Panels are labeled by the real and imaginary parts of the eigenvalues. 
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to provide much more detail, essentially because we are dealing with continuous discs rather 
than a finite number of rings. 

As 77 incre ases from 0, t he eigenvalues go from real to complex, a bifurcation that has 



been traced in iToumal ( 12002| ) to a phenomenon identified by M. J. Krein due to the resonant 
crossing of stable modes. The complex eigenvalues come in complex conjugate pairs, so there 
are two branches to the distribution. As rj increases, the branches progressively separate and, 
for ?7 = 1/2 must lie along the positive and negative imaginary axes. It is intriguing that 
each of these two branches consists of more than one arm. In the close-up views provided 
by the (b) panels, it appears as if each of the branches has two arms; however, more detailed 
investigations are required to determine if there are more arms. The arms of each of the 
branches are most widely separated when rj = 0.25, which is the value of rj exactly midway 
in its range ^ 77 ^ 0.5 . The separations decrease as rj approaches either or 1/2; this is 
natural because, for rj = both branches must lie on the real axis and, for rj = 1/2 both 
branches must lie on the imaginary axis. 

The eigenfunctions are in general complex, and have a rich structure as functions of their 
eigenvalues. Since our interest is in the unstable modes, we now display in Fig. [7| plots of 
the Zji = ^[Z{R)] in Eq. fl29l) . corresponding to the "most unstable" modes (for softening 
parameter b = 0.1 and for t] equal to 0.1, 0.25 and 0.4). In other words, for some chosen 
value of aR, we display the real part of the eigenf unction corresponding to the largest value 
of aj. For a fixed value of 77, the number of nodes of the eigenfunctions decreases with 
increasing pattern speed and growth rate. Fig [8] is a gray-scale plot of the surface density 
perturbations in the ± discs, S^(i?, 0, t = 0), for the parameter values rj = 0.25 and 6 = 0.1, 
and a = 0.1408 + 10.2332 . 

It is also of interest to ask how eigenfunctions from two different arms of the same branch 
behave. To do this, we picked two eigenfunctions with nearly the same value of a^, but with 
values of aj corresponding to two different arms of one branch; Fig. [9] shows two such pairs of 
eigenfunctions for 77 = 0.25. We have looked at pairs of such eigenfunctions for other values 
of ?7, but do not display them, here we note what seems to be a general trend: (i) the two 
members of a pair are more similar to each other when the values of their an are closer to 
each other; (ii) the member of a pair with the smaller value of ai is more displaced toward 
larger radii. 
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5 CONCLUSIONS 



We study linear, slow m = 1 modes in softened gravity, counter-rotating Keplerian discs. 
The eigenvalue problem is formulated as a pair of coupled integral equations for the ± 
modes. We then specialize to the case when the two discs have similar surface density pro- 
files but different ± disc masses. It is of great interest to study the properties of the modes 
as a function of rj, which is the fraction of the total disc mass in the retrograde population. 
Recasting the coupled equations as a single equation in a new modal variable, we are able 
to demonstrate some general properties: for instance, when rj = 1/2, the eigenvalues must 
be purely imaginary or, equivalently, the modes are purely unstable. In other words, when 
the ± discs have identical surface density profiles then t here are grow i ng m = 1 modes with 
zero p a ttern speed, a conclu s ion w h ich is consistent with 



fll990h : 



Sellwood fc Merritt 



Araki](ll987f): 



fll994[ ): iLovelace et all (119971 1: iToumal ( 120021 ): iTremaind (12005f ). 



Palmer fc Papaloizou 



To study modes for general values of rj, the eigenva l ue pr oblem needs to be solved numer- 
ically. Our method is broadly based on iTremaind ( 120011 ). but there are some differences 
whose details have been discussed i n the text . The main point of departure is in the way 



that softening has been treated. In 



Tremaind (1200 ll ). a dimensionless softening parameter 



/3 = b/R was introduced, and the eigenvalue problem was solved by holding the parameter 
P constant. This procedure renders the physical softening length, b, effectively dependent on 
radius (making it larger at larger radii), thereby not corresponding to any simple force law 
between two disc particles located at different radii. We have preferred to keep b constant, 
so that the force law between two disc particles is through the usual Miller prescription. 

We calculate eigenvalues and eigenf unctions numerically for discs with surface density 
profiles of Kuzmin form. Kuzmin discs, being centrally concentrated, are r easonable candi 



dates for unperturbed discs. Moreover, earlier investigations of slow modes (JTremaine 



Sridhar fc Saini 



2001 : 



2010l ) have explored modes in Kuzmin discs, so this choice is p articularly 



usefu 



Tremaine 



for comparisons with earlier work. Comparing our results with those of 
( 1200 ll ) for 1] = (when the slow modes are stable), we find that the eigenfunctions are of 
broadly similar form, but the eigenvalues differ by up to ~ 30% ; this is a result of the differ- 
ent ways in which we have treated softening. For the case of no net rotation (77 = 1/2), we find 
that the growth rates (of the unstable modes) we calculate are broadly co nsistent with the 



Sridhar fc Saini 



short- wavelength branch of the global WKB modes determined earlier by 

(I2OIOI ). but not their long-wavelength branch. This disagreement probably arises because 
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the long- wavelength branch corresponds to wavelengths of order the disc scale length where 
WKB approximation breaks down. Moreover, the agreement between our results and their 
short-wavelength branch holds only in a broad sense, because there are differences in the 
numerical values of the eigenvalues. We trace this difference to the fact that 



Sridhar fc Saini 



( I2OIOI ) used an analytical result for the precession frequency corresponding to unsoftened 
gravity, whereas we have consistently used softened gravity for all gravitational interactions 
between disc particles. 

We have also investigate eigenmodes for values of t] other than and 1/2. These cases ar e 



Sridhar &: Saini 



mm, 



particularly interesting, not only because they were not explored by 
but because the eigenvalues can be truly complex, corresponding to growing (and damped) 
modes with non zero pattern speeds. We have presented results for rj = 0.1,0.25 and 0.4 
in the previous sections. Based on these, we interpolate and offer the following conclusions 
about the properties of the eigenmodes and their physical implications, for all values of r) 
(which is the mass fraction in the retrograde population): 

(i) For a general value of 77 (between and 1/2), the distribution of eigenvalues in the 
complex plane has two branches. These branches are symmetrically placed about the real 
axis, because the eigenvalues come in complex conjugate pairs. 

(ii) The pattern speed appears to be non negative for all values of rj, with the growth (or 
damping) rate being larger for larger values of the pattern speed. 

(iii) For a fixed value of rj, the number of nodes of the eigenfunctions decreases with 
increasing pattern speed and growth (or damping) rate. 

(iv) For a value of pattern speed in a chosen narrow interval, the growth (or damping) 
rate increases as rj increases from to 1/2. 

(v) Each of the two branches in the complex (eigenvalue) plane has at least two arms. 
When rj = 0, the eigenvalues are all real, so both branches lie on the real axis, with zero 
spacing between the arms. As 77 increases, the branches lift out of the real axis, and the 
arms separate. It appears as if the maximum separation between the arms happens when 
rj = 1/4. As rj increases further, the branches continue to rise with greater slope, while the 
arm separation begins decreasing. Finally, when 77 = 1/2, the arm separation decreases to 
zero as the branches lie on the imaginary axis. 

Observations of lopsided brightness distributions around massive black holes are some- 
what more likely to favour the detection of modes with fewer nodes than modes with a large 
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number of nodes, because the former suffer less cancellation due to finite angular resolution. 
From items (ii) and (iii) above, we note that the modes with a small number of nodes also 
happen to be those with larger values of the pattern speed and growth rate, both qual- 
ities that enable detection to a greater degree. Having said this, it would be appropriate 
to note some limitations of our work. Softened gravity discs are, after all, surrogates for 
discs composed of collisionless particles (such as stars) with non zero thickness and velocity 
dispersions. It is necessary to formulate the eigenvalue problem for truly collisionless discs, 
in order to really deal with stellar discs around massive black holes. Meanwhile, our results 
will serve as a benchmark for future investigations of modes in these more realistic models. 
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APPENDIX A: EXPRESSING SOFTENED LAPLACE COEFFICIENTS IN 
TERMS OF (UNSOFTENED) LAPLACE COEFFICIENTS 

Softened Laplace coefficients were defined in 



Touma 



(120021 ) as 



where 



We now write 



One solution for 7 and 6 is. 



2 r 

B^{a,P) = - de 

TT Jo 



COS m9 



A = l + a^ + P^ -2acos9 



A = 7" + (5' -275 cos ^ 



(Al) 



(A2) 



(A3) 



7 



l + a^ + Z^^ , 1 /-— — „ , ^„,^ -—, 

+ - V(l + a^ + /32)2 - 4a2 



1/2 



6 = ^ 

7 

Therefore 



(A4) 



i?r(«,/3) = l-'bT/2iSh) 



(A5) 
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where 



&r/2(«) 



- I de- 



ll 



cos mO 



:i + «2 



2a cos 0) 



s/2 ' 



a < 1 



(A6) 



are the famihar (unsoftened) Laplace coefficients (jMurray &: Derniottlll999l ). From eqn. (1A5I1 . 
we must have ((5/7) < 1. That this is indeed true can be proved using eqns. (JA4]): 7 is a 
monotonically increasing function of /3^, hence 7^1, and (5/7) = (a/7^) < 1 . 
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